module obs_operator_interp_mod

  use kdtree
  use synop_mod
  use raob_mod
  use amdar_mod
  use container
  use const_mod
  use sphere_geometry_mod
  use efso_utils_mod
  use obs_synop_mod
  use obs_raob_mod
  use obs_amdar_mod
  use obs_profiler_mod
  use latlon_mesh_mod
  use latlon_data_mod
  use latlon_parallel_mod
  use latlon_process_mod
  use obs_operator_types_mod

  implicit none

  private

  public obs_operator_interp_prepare
  public obs_operator_interp_run

  interface obs_operator_interp_prepare
    module procedure obs_operator_interp_prepare_latlon_mesh
  end interface obs_operator_interp_prepare

  interface obs_operator_interp_run
    module procedure obs_operator_interp_run_latlon_mesh
  end interface obs_operator_interp_run

contains

  subroutine obs_operator_interp_prepare_latlon_mesh(obs, mesh)

    class(*), intent(inout), target :: obs
    type(latlon_mesh_type), intent(in) :: mesh

    type(kdtree_type) tree
    type(linked_list_iterator_type) it
    type(obs_grid_connector_type), pointer :: c
    integer, allocatable :: ig(:)
    integer, allocatable :: jg(:)
    real(r8), allocatable :: xg(:,:)
    real(r8) x(3)
    integer i, j, l

    ! 1. Create KDTree for mesh grids.
    allocate(xg(3,mesh%num_lon*mesh%num_lat))
    allocate(ig(mesh%num_lon*mesh%num_lat))
    allocate(jg(mesh%num_lon*mesh%num_lat))
    l = 1
    do j = mesh%lat_ibeg, mesh%lat_iend
      do i = mesh%lon_ibeg, mesh%lon_iend
        call cartesian_transform(mesh%lon(i), mesh%lat(j), x(1), x(2), x(3))
        xg(:,l) = x; ig(l) = i; jg(l) = j
        l = l + 1
      end do
    end do
    call tree%build(xg)

    ! 2. Find the enclosing grids for each obs point.
    select type (obs)
    type is (obs_synop_type)
      it = linked_list_iterator(obs%records); i = 1
      do while (.not. it%ended())
        c => obs%connectors(i)
        select type (record => it%value)
        type is (synop_record_type)
          call c%init(4, 2, record)
          call find_neighbors(mesh, tree, ig, jg,  record%station%lon * rad, record%station%lat * rad, c)
        end select
        call it%next(); i = i + 1
      end do
    type is (obs_raob_type)
      it = linked_list_iterator(obs%records); i = 1
      do while (.not. it%ended())
        select type (record => it%value)
        type is (raob_record_type)
          ! Mandatory levels
          c => obs%connectors_man(i)
          call c%init(4, 2, record, obs_plev=record%man%p)
          call find_neighbors(mesh, tree, ig, jg,  record%station%lon * rad, record%station%lat * rad, c)
          ! Significant temperature levels
          c => obs%connectors_sigt(i)
          call c%init(4, 2, record, obs_plev=record%sigt%p)
          call find_neighbors(mesh, tree, ig, jg,  record%station%lon * rad, record%station%lat * rad, c)
          ! Significant wind levels
          c => obs%connectors_sigw(i)
          call c%init(4, 2, record, obs_plev=record%sigw%p)
          call find_neighbors(mesh, tree, ig, jg,  record%station%lon * rad, record%station%lat * rad, c)
          ! Tropopause levels
          c => obs%connectors_trop(i)
          call c%init(4, 2, record, obs_plev=record%trop%p)
          call find_neighbors(mesh, tree, ig, jg,  record%station%lon * rad, record%station%lat * rad, c)
        end select
        call it%next(); i = i + 1
      end do
    type is (obs_amdar_type)
      it = linked_list_iterator(obs%records); i = 1
      do while (.not. it%ended())
        c => obs%connectors(i)
        select type (record => it%value)
        type is (amdar_record_type)
          call c%init(4, 2, record, obs_alt=record%h, obs_plev=[record%p])
          call find_neighbors(mesh, tree, ig, jg, record%lon * rad, record%lat * rad, c)
        end select
        call it%next(); i = i + 1
      end do
    type is (obs_profiler_type)
      it = linked_list_iterator(obs%records); i = 1
      do while (.not. it%ended())
        select type (record => it%value)
        type is (profiler_record_type)
          c => obs%connectors(i)
          call c%init(4, 2, record, obs_zlev=record%pro%h)
          call find_neighbors(mesh, tree, ig, jg,  record%station%lon * rad, record%station%lat * rad, c)
        end select
        call it%next(); i = i + 1
      end do
    end select

    deallocate(xg, ig, jg)

  end subroutine obs_operator_interp_prepare_latlon_mesh

  subroutine obs_operator_interp_run_latlon_mesh(mesh, model_data, iens, obs)

    type(latlon_mesh_type), intent(in) :: mesh
    type(latlon_data_type), intent(in) :: model_data
    integer, intent(in) :: iens
    class(*), intent(inout), target :: obs

    type(obs_grid_connector_type), pointer :: c
    real(r8) vert_wgt(2), zm(mesh%lev_lb:mesh%lev_ub), dz, dlnp
    integer i, imask, k, vert_idx(2)

    imask = 1
    select type (obs)
    type is (obs_synop_type)
      do i = 1, size(obs%connectors)
        c => obs%connectors(i)
        if (iens == -1) then
          obs%simulated_mean_records(i)%ua = bilinear_2d(mesh, model_data%u10, c%grid_idx, c%grid_wgt)
          obs%simulated_mean_records(i)%va = bilinear_2d(mesh, model_data%v10, c%grid_idx, c%grid_wgt)
          obs%simulated_mean_records(i)%ta = bilinear_2d(mesh, model_data%t2 , c%grid_idx, c%grid_wgt)
          obs%simulated_mean_records(i)%p  = bilinear_2d(mesh, model_data%ps , c%grid_idx, c%grid_wgt)
          obs%simulated_mean_records(i)%sh = bilinear_2d(mesh, model_data%q2 , c%grid_idx, c%grid_wgt)
        else
          obs%simulated_records(i,iens)%ua = bilinear_2d(mesh, model_data%u10, c%grid_idx, c%grid_wgt)
          obs%simulated_records(i,iens)%va = bilinear_2d(mesh, model_data%v10, c%grid_idx, c%grid_wgt)
          obs%simulated_records(i,iens)%ta = bilinear_2d(mesh, model_data%t2 , c%grid_idx, c%grid_wgt)
          obs%simulated_records(i,iens)%p  = bilinear_2d(mesh, model_data%ps , c%grid_idx, c%grid_wgt)
          obs%simulated_records(i,iens)%sh = bilinear_2d(mesh, model_data%q2 , c%grid_idx, c%grid_wgt)
        end if
        ! Correct to real terrain height.
        ! ...
      end do
    type is (obs_raob_type)
      do i = 1, obs%records%size
        if (iens == -1) then
          call obs%mean_simulated_records(i)%init_with_size(num_man_level =size(obs%connectors_man (i)%obs_plev), &
                                                            num_sigt_level=size(obs%connectors_sigt(i)%obs_plev), &
                                                            num_sigw_level=size(obs%connectors_sigw(i)%obs_plev), &
                                                            num_trop_level=size(obs%connectors_trop(i)%obs_plev))
          call obs%simulated_mean_records(i)%init_with_size(num_man_level =size(obs%connectors_man (i)%obs_plev), &
                                                            num_sigt_level=size(obs%connectors_sigt(i)%obs_plev), &
                                                            num_sigw_level=size(obs%connectors_sigw(i)%obs_plev), &
                                                            num_trop_level=size(obs%connectors_trop(i)%obs_plev))
        else
          call obs%simulated_records(i,iens)%init_with_size(num_man_level =size(obs%connectors_man (i)%obs_plev), &
                                                            num_sigt_level=size(obs%connectors_sigt(i)%obs_plev), &
                                                            num_sigw_level=size(obs%connectors_sigw(i)%obs_plev), &
                                                            num_trop_level=size(obs%connectors_trop(i)%obs_plev))
        end if
        ! Mandatory levels
        c => obs%connectors_man(i)
        do k = 1, size(c%obs_plev)
          ! Find out the model pressure level that equals the obs pressure level.
          vert_idx = ind(mesh%lev, c%obs_plev(k))
          if (vert_idx(1) == -999 .and. vert_idx(2) == -999) then
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
          else
            vert_wgt(1) = (mesh%lev(vert_idx(2)) - c%obs_plev(k)) / &
                          (mesh%lev(vert_idx(2)) - mesh%lev(vert_idx(1)))
            vert_wgt(2) = 1 - vert_wgt(1)
            if (iens == -1) then
              obs%simulated_mean_records(i)%man%p (k) = c%obs_plev(k)
              obs%simulated_mean_records(i)%man%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%man%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%man%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%man%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            else
              obs%simulated_records(i,iens)%man%p (k) = c%obs_plev(k)
              obs%simulated_records(i,iens)%man%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%man%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%man%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%man%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              !select type (r => c%record)
              !type is (raob_record_type)
              !  print *, i, iens, k, &
              !    obs%simulated_records(i,iens)%man%ua(k), r%man%ua(k), &
              !    obs%simulated_records(i,iens)%man%va(k), r%man%va(k), &
              !    obs%simulated_records(i,iens)%man%ta(k), r%man%ta(k), &
              !    obs%simulated_records(i,iens)%man%p (k), r%man%p (k), &
              !    obs%simulated_records(i,iens)%man%sh(k), r%man%sh(k)
              !end select
            end if
          end if
        end do
        ! Significant temperature levels
        c => obs%connectors_sigt(i)
        do k = 1, size(c%obs_plev)
          ! Find out the model pressure level that equals the obs pressure level.
          vert_idx = ind(mesh%lev, c%obs_plev(k))
          if (vert_idx(1) == -999 .or. vert_idx(2) == -999) then
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
          else
            vert_wgt(1) = (mesh%lev(vert_idx(2)) - c%obs_plev(k)) / &
                          (mesh%lev(vert_idx(2)) - mesh%lev(vert_idx(1)))
            vert_wgt(2) = 1 - vert_wgt(1)
            if (iens == -1) then
              obs%simulated_mean_records(i)%sigt%p (k) = c%obs_plev(k)
              obs%simulated_mean_records(i)%sigt%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%sigt%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%sigt%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%sigt%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            else
              obs%simulated_records(i,iens)%sigt%p (k) = c%obs_plev(k)
              obs%simulated_records(i,iens)%sigt%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%sigt%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%sigt%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%sigt%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              !select type (r => c%record)
              !type is (raob_record_type)
              !  print *, i, iens, k, &
              !    obs%simulated_records(i,iens)%sigt%ua(k), r%sigt%ua(k), &
              !    obs%simulated_records(i,iens)%sigt%va(k), r%sigt%va(k), &
              !    obs%simulated_records(i,iens)%sigt%ta(k), r%sigt%ta(k), &
              !    obs%simulated_records(i,iens)%sigt%p (k), r%sigt%p (k), &
              !    obs%simulated_records(i,iens)%sigt%sh(k), r%sigt%sh(k)
              !end select
            end if
          end if
        end do
        ! Significant wind levels
        c => obs%connectors_sigw(i)
        do k = 1, size(c%obs_plev)
          ! Find out the model pressure level that equals the obs pressure level.
          vert_idx = ind(mesh%lev, c%obs_plev(k))
          if (vert_idx(1) == -999 .or. vert_idx(2) == -999) then
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
          else
            vert_wgt(1) = (mesh%lev(vert_idx(2)) - c%obs_plev(k)) / &
                          (mesh%lev(vert_idx(2)) - mesh%lev(vert_idx(1)))
            vert_wgt(2) = 1 - vert_wgt(1)
            if (iens == -1) then
              obs%simulated_mean_records(i)%sigw%p (k) = c%obs_plev(k)
              obs%simulated_mean_records(i)%sigw%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%sigw%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%sigw%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%sigw%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            else
              obs%simulated_records(i,iens)%sigw%p (k) = c%obs_plev(k)
              obs%simulated_records(i,iens)%sigw%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%sigw%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%sigw%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%sigw%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              !select type (r => c%record)
              !type is (raob_record_type)
              !  print *, i, iens, k, &
              !    obs%simulated_records(i,iens)%sigw%ua(k), r%sigw%ua(k), &
              !    obs%simulated_records(i,iens)%sigw%va(k), r%sigw%va(k), &
              !    obs%simulated_records(i,iens)%sigw%ta(k), r%sigw%ta(k), &
              !    obs%simulated_records(i,iens)%sigw%p (k), r%sigw%p (k), &
              !    obs%simulated_records(i,iens)%sigw%sh(k), r%sigw%sh(k)
              !end select
            end if
          end if
        end do
        ! Tropopause levels
        c => obs%connectors_trop(i)
        do k = 1, size(c%obs_plev)
          ! Find out the model pressure level that equals the obs pressure level.
          vert_idx = ind(mesh%lev, c%obs_plev(k))
          if (vert_idx(1) == -999 .or. vert_idx(2) == -999) then
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
            obs%mask(imask) = 0; imask = imask + 1
          else
            vert_wgt(1) = (mesh%lev(vert_idx(2)) - c%obs_plev(k)) / &
                          (mesh%lev(vert_idx(2)) - mesh%lev(vert_idx(1)))
            vert_wgt(2) = 1 - vert_wgt(1)
            if (iens == -1) then
              obs%simulated_mean_records(i)%trop%p (k) = c%obs_plev(k)
              obs%simulated_mean_records(i)%trop%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%trop%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%trop%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_mean_records(i)%trop%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            else
              obs%simulated_records(i,iens)%trop%p (k) = c%obs_plev(k)
              obs%simulated_records(i,iens)%trop%ua(k) = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%trop%va(k) = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%trop%ta(k) = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              obs%simulated_records(i,iens)%trop%sh(k) = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
              !select type (r => c%record)
              !type is (raob_record_type)
              !  print *, i, iens, k, &
              !    obs%simulated_records(i,iens)%trop%ua(k), r%trop%ua(k), &
              !    obs%simulated_records(i,iens)%trop%va(k), r%trop%va(k), &
              !    obs%simulated_records(i,iens)%trop%ta(k), r%trop%ta(k), &
              !    obs%simulated_records(i,iens)%trop%p (k), r%trop%p (k), &
              !    obs%simulated_records(i,iens)%trop%sh(k), r%trop%sh(k)
              !end select
            end if
          end if
        end do
      end do
      call global_min(proc%comm, obs%mask)
    type is (obs_amdar_type)
      do i = 1, size(obs%connectors)
        c => obs%connectors(i)
        ! Find out the model pressure level that equals the obs pressure level.
        vert_idx = ind(mesh%lev, c%obs_plev(1))
        if (vert_idx(1) == -999 .or. vert_idx(2) == -999) then
          obs%mask(imask) = 0; imask = imask + 1
          obs%mask(imask) = 0; imask = imask + 1
          obs%mask(imask) = 0; imask = imask + 1
          obs%mask(imask) = 0; imask = imask + 1
        else
          vert_wgt(1) = (mesh%lev(vert_idx(2)) - c%obs_plev(1)) / &
                        (mesh%lev(vert_idx(2)) - mesh%lev(vert_idx(1)))
          vert_wgt(2) = 1 - vert_wgt(1)
          if (iens == -1) then
            obs%simulated_mean_records(i)%z  = c%obs_alt
            obs%simulated_mean_records(i)%p  = c%obs_plev(1)
            obs%simulated_mean_records(i)%ua = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            obs%simulated_mean_records(i)%va = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            obs%simulated_mean_records(i)%ta = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            obs%simulated_mean_records(i)%sh = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
          else
            obs%simulated_records(i,iens)%z  = c%obs_alt
            obs%simulated_records(i,iens)%p  = c%obs_plev(1)
            obs%simulated_records(i,iens)%ua = bilinear_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            obs%simulated_records(i,iens)%va = bilinear_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            obs%simulated_records(i,iens)%ta = bilinear_3d(mesh, model_data%t, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
            obs%simulated_records(i,iens)%sh = bilinear_3d(mesh, model_data%q, c%grid_idx, c%grid_wgt, vert_idx, vert_wgt); imask = imask + 1
          end if
        end if
      end do
      call global_min(proc%comm, obs%mask)
    type is (obs_profiler_type)
      if (iens == -1 .and. proc%is_root_proc()) call log_notice('Interpolate profiler pressure using ensemble mean data.')
      do i = 1, size(obs%connectors)
        if (iens == -1) then
          call obs%mean_simulated_records(i)%init_with_size(num_level=size(obs%connectors(i)%obs_zlev))
          call obs%simulated_mean_records(i)%init_with_size(num_level=size(obs%connectors(i)%obs_zlev))
        else
          call obs%simulated_records(i,iens)%init_with_size(num_level=size(obs%connectors(i)%obs_zlev))
        end if
        c => obs%connectors(i)
        call bilinear_profile(mesh, model_data%z, c%grid_idx, c%grid_wgt, zm)
        do k = 1, size(c%obs_zlev)
          ! Find out the model pressure level that equals the obs pressure level.
          if (iens == -1) then
            obs%simulated_mean_records(i)%pro%h (k) = c%obs_zlev(k)
            obs%simulated_mean_records(i)%pro%ua(k) = bilinear_z_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, model_data%z, c%obs_zlev(k))
            if (is_missing(obs%simulated_mean_records(i)%pro%ua(k))) obs%mask(imask) = 0; imask = imask + 1
            obs%simulated_mean_records(i)%pro%va(k) = bilinear_z_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, model_data%z, c%obs_zlev(k))
            if (is_missing(obs%simulated_mean_records(i)%pro%va(k))) obs%mask(imask) = 0; imask = imask + 1
            ! Interpolate profiler pressure using model.
            vert_idx = ind(zm, c%obs_zlev(k))
            if (vert_idx(1) /= -999 .and. vert_idx(2) /= -999) then
              select type (r => c%record)
              type is (profiler_record_type)
                dz = zm(vert_idx(2)) - zm(vert_idx(1))
                dlnp = log(mesh%lev(vert_idx(2))) - log(mesh%lev(vert_idx(1)))
                r%pro%p(k) = exp(log(mesh%lev(vert_idx(1))) + (dlnp / dz) * (c%obs_zlev(k) - zm(vert_idx(1))))
                !if (proc%is_root_proc()) then
                !  print *, i, k, r%pro%h(k), r%pro%p(k), mesh%lev(vert_idx(1)), dz, dlnp, c%obs_zlev(k), zm(vert_idx(1))
                !end if
              end select
            end if
          else
            obs%simulated_records(i,iens)%pro%h (k) = c%obs_zlev(k)
            obs%simulated_records(i,iens)%pro%ua(k) = bilinear_z_3d(mesh, model_data%u, c%grid_idx, c%grid_wgt, model_data%z, c%obs_zlev(k))
            if (is_missing(obs%simulated_records(i,iens)%pro%ua(k))) obs%mask(imask) = 0; imask = imask + 1
            obs%simulated_records(i,iens)%pro%va(k) = bilinear_z_3d(mesh, model_data%v, c%grid_idx, c%grid_wgt, model_data%z, c%obs_zlev(k))
            if (is_missing(obs%simulated_records(i,iens)%pro%va(k))) obs%mask(imask) = 0; imask = imask + 1
            !if (proc%id == 26 .and. i == 74) then
            !  select type (r => c%record)
            !  type is (profiler_record_type)
            !    print *, trim(r%station%name), c%obs_zlev(k), &
            !      obs%simulated_records(i,iens)%pro%ua(k), r%pro%ua(k), obs%mask(imask-2), &
            !      obs%simulated_records(i,iens)%pro%va(k), r%pro%va(k), obs%mask(imask-1)
            !  end select
            !end if
          end if
        end do
      end do
      call global_min(proc%comm, obs%mask)
    end select

  end subroutine obs_operator_interp_run_latlon_mesh

  subroutine find_neighbors(mesh, tree, ig, jg,  lon, lat, connector)

    type(latlon_mesh_type), intent(in) :: mesh
    type(kdtree_type), intent(in) :: tree
    integer, intent(in) :: ig(:)
    integer, intent(in) :: jg(:)
    real(r8), intent(in) :: lon
    real(r8), intent(in) :: lat
    type(obs_grid_connector_type), intent(inout) :: connector

    real(r8) x(3), a, b
    integer ngb_idx(1), i0, j0

    call cartesian_transform(lon, lat, x(1), x(2), x(3))
    call tree%search(x, ngb_idx)
    connector%obs_lon = lon
    connector%obs_lat = lat
    i0 = ig(ngb_idx(1))
    j0 = jg(ngb_idx(1))
    if (lon >= mesh%lon(i0)) then
      connector%grid_idx(1,1) = i0
      connector%grid_idx(1,2) = i0 + 1
      connector%grid_idx(1,3) = i0 + 1
      connector%grid_idx(1,4) = i0
    else
      connector%grid_idx(1,1) = i0 - 1
      connector%grid_idx(1,2) = i0
      connector%grid_idx(1,3) = i0
      connector%grid_idx(1,4) = i0 - 1
    end if
    if (lat >= mesh%lat(j0)) then
      connector%grid_idx(2,1) = j0
      connector%grid_idx(2,2) = j0
      connector%grid_idx(2,3) = j0 + 1
      connector%grid_idx(2,4) = j0 + 1
    else
      connector%grid_idx(2,1) = j0 - 1
      connector%grid_idx(2,2) = j0 - 1
      connector%grid_idx(2,3) = j0
      connector%grid_idx(2,4) = j0
    end if
    if (connector%grid_idx(1,1) < 1) then
      connector%around_zero_lon = .true.
      connector%grid_idx(1,1) = connector%grid_idx(1,1) + mesh%num_lon + 1
      connector%grid_idx(1,4) = connector%grid_idx(1,4) + mesh%num_lon + 1
    else if (connector%grid_idx(1,2) > mesh%num_lon) then
      connector%around_zero_lon = .true.
      connector%grid_idx(1,2) = connector%grid_idx(1,2) - mesh%num_lon
      connector%grid_idx(1,3) = connector%grid_idx(1,3) - mesh%num_lon
    end if

    if (connector%around_zero_lon .and. lon < pi2 - lon) then
      ! Western Sphere
      a = (lon - (mesh%lon(connector%grid_idx(1,1)) - pi2)) / &
          (mesh%lon(connector%grid_idx(1,2)) - mesh%lon(connector%grid_idx(1,1)) + pi2)
    else
      a = (lon - mesh%lon(connector%grid_idx(1,1))) / (mesh%lon(connector%grid_idx(1,2)) - mesh%lon(connector%grid_idx(1,1)))
    end if
    b = (lat - mesh%lat(connector%grid_idx(2,1))) / (mesh%lat(connector%grid_idx(2,4)) - mesh%lat(connector%grid_idx(2,1)))

    connector%grid_wgt(1) = (1 - a) * (1 - b)
    connector%grid_wgt(2) =      a  * (1 - b)
    connector%grid_wgt(3) =      a  *      b
    connector%grid_wgt(4) = (1 - a) *      b

  end subroutine find_neighbors

  real(r8) function bilinear_2d(mesh, x, idx, wgt) result(res)

    type(latlon_mesh_type), intent(in) :: mesh
    real(r8), intent(in) :: x(mesh%lon_lb:mesh%lon_ub,mesh%lat_lb:mesh%lat_ub)
    integer, intent(in) :: idx(2,4)
    real(r8), intent(in) :: wgt(4)

    res = x(idx(1,1),idx(2,1)) * wgt(1) + &
          x(idx(1,2),idx(2,2)) * wgt(2) + &
          x(idx(1,3),idx(2,3)) * wgt(3) + &
          x(idx(1,4),idx(2,4)) * wgt(4)

  end function bilinear_2d

  real(r8) function bilinear_3d(mesh, x, hor_idx, hor_wgt, vert_idx, vert_wgt) result(res)

    type(latlon_mesh_type), intent(in) :: mesh
    real(r8), intent(in) :: x(mesh%lon_lb:mesh%lon_ub,mesh%lat_lb:mesh%lat_ub,mesh%lev_lb:mesh%lev_ub)
    integer, intent(in) :: hor_idx(2,4)
    real(r8), intent(in) :: hor_wgt(4)
    integer, intent(in) :: vert_idx(2)
    real(r8), intent(in) :: vert_wgt(2)

    real(r8) tmp(4)
    integer i

    do i = 1, 4
      tmp(i) = vert_wgt(1) * x(hor_idx(1,i),hor_idx(2,i),vert_idx(1)) + &
               vert_wgt(2) * x(hor_idx(1,i),hor_idx(2,i),vert_idx(2))
    end do
    res = sum(tmp * hor_wgt)

  end function bilinear_3d

  real(r8) function bilinear_z_3d(mesh, x, hor_idx, hor_wgt, grid_z, z) result(res)

    type(latlon_mesh_type), intent(in) :: mesh
    real(r8), intent(in) :: x(mesh%lon_lb:mesh%lon_ub,mesh%lat_lb:mesh%lat_ub,mesh%lev_lb:mesh%lev_ub)
    integer, intent(in) :: hor_idx(2,4)
    real(r8), intent(in) :: hor_wgt(4)
    real(r8), intent(in) :: grid_z(mesh%lon_lb:mesh%lon_ub,mesh%lat_lb:mesh%lat_ub,mesh%lev_lb:mesh%lev_ub)
    real(r8), intent(in) :: z

    real(r8) z1, z2, tmp(4)
    integer i, vert_idx(2)

    do i = 1, 4
      vert_idx = ind(grid_z(hor_idx(1,i),hor_idx(2,i),:), z)
      if (vert_idx(1) == -999 .or. vert_idx(2) == -999) then
        res = real_missing_value
        return
      else
        z1 = grid_z(hor_idx(1,i),hor_idx(2,i),vert_idx(1))
        z2 = grid_z(hor_idx(1,i),hor_idx(2,i),vert_idx(2))
        tmp(i) = ((z2 - z) * x(hor_idx(1,i),hor_idx(2,i),vert_idx(1)) + &
                  (z - z1) * x(hor_idx(1,i),hor_idx(2,i),vert_idx(2))) / (z2 - z1)
      end if
    end do
    res = sum(tmp * hor_wgt)

  end function bilinear_z_3d

  subroutine bilinear_profile(mesh, x, hor_idx, hor_wgt, y)

    type(latlon_mesh_type), intent(in) :: mesh
    real(r8), intent(in) :: x(mesh%lon_lb:mesh%lon_ub,mesh%lat_lb:mesh%lat_ub,mesh%lev_lb:mesh%lev_ub)
    integer, intent(in) :: hor_idx(2,4)
    real(r8), intent(in) :: hor_wgt(4)
    real(r8), intent(out) :: y(mesh%lev_lb:mesh%lev_ub)

    real(r8) tmp(4)
    integer i, k, vert_idx(2)

    do k = mesh%lev_lb, mesh%lev_ub
      do i = 1, 4
        tmp(i) = x(hor_idx(1,i),hor_idx(2,i),k)
      end do
      y(k) = sum(tmp * hor_wgt)
    end do

  end subroutine bilinear_profile

end module obs_operator_interp_mod
